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SUMMARY 

The  analysis  of  random  data  is  frequently  carried  out  'on-line’  using  a  hardware 
based  Fourier  analyser.  Capabilities  similar  to  those  provided  in  the  ‘on-line’  analysers 
have  been  developed  in  a  Fortran  computer  program  that  can  be  used  to  analyse  a  digital 
time-dependent  sigrud  at  the  post-test  stage.  An  overview  of  the  calculation  methods  used 
within  this  software  is  presented. 
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Notation 


c  Modified  boxcar  fimctiou 

/  frequency 

Sc  NyqoUt  frequency 

G  Power  spectral  density 

E  Transfer  fnnction 

h  Sampling  time  interval 

I  Number  of  a4jacent  frequencies  used  in  frequency  smoothing 

N  Number  of  discrete  points  used  in  each  FFT  computation 

q  Number  of  ensembles  used  in  the  averapng  process 
R  Correlation  function 

s/  Scaie  factor 

T  Signal  length 

t  Time 

U  Fourier  transform  of  u 

u  Boxcar  function 

X,  Y  Fourier  transforms  of  x  and  y 
X,  y  Input  signsds 

7  Coherence  function 

6  Argument  of  cross  power  spectral  density  function 

r  Time  lag 

^  Phase  of  transfer  function 

9  Mean  square  value 

Superscripts 

~  Raw  estimate 
~  Smoothed  estimate 
'  Modified  function 
*  Compiex  conjugate 
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1  Introduction 


The  deTelopment  of  the  Fast  Fourier  Traniform  algorithm  has  extended  the  procedures 
available  for  the  analysis  of  random  data.  These  procedures  can  be  carried  out  on-line 
using  a  hardware  based  machine,  or  the  data  can  be  processed  using  computer  software. 

The  Wavetek  804A  machine  is  an  example  of  a  digital  Fourier  analysis  system  which 
can  process  unsteady  signab  ‘real  time’.  To  provide  post-test  analysis  facilities  with 
similar  capabilities,  a  PC-based  computer  program  has  been  developed.  The  computer 
program  SPEC  which  has  been  written  in  Fortran  77  uses  a  Fast  Fourier  lyansfarm 
routine  as  the  basis  for  aU  subsequent  calculations  and  can  be  linked  to  the  public 
domain  program  GNUPLOT  for  graphical  output. 

This  report  introduces  the  user  to  basic  spectral  analysis  theory  and  highlights  some 
aspects  of  data  sampling  and  interpretation  that  need  to  be  considered  when  dealing 
with  discrete  Fourier  transforms.  A  more  detailed  description  of  spectral  analysis  is 
given  by  Bendat  and  Piersol  (Reference  [1]).  This  text  forms  the  basis  for  much  of  the 
theory  presented  in  this  report. 

An  outline  of  the  computer  program  SPEC  is  also  given  with  descriptions  of  both 
input  and  output. 

2  The  Fourier  Integral  and  Discrete  Fourier  Transform 

In  this  report,  the  equations  used  in  the  analysis  program  are  presented.  The  basic 
properties  and  classification  of  random  data,  its  measurement  and  analysis  are  described 
in  Reference  [1].  Consideration  of  the  data  characteristics  should  be  made  prior  to  the 
selection  of  the  analysis  techniques. 

In  the  analysis  of  transient  data,  the  Fourier  Integral  is  lued  to  determine  a  frequency 
spectrum  firom  any  arbitrary  time  dependent  signal.  The  infinite  range  Fourier  Integral 
is  given  by  (Ref.  [1]): 

-*■(/)  =  r  (1) 

V— oo 

Realistically,  the  signal  of  interest  is  defined  over  a  finite  time  interval  [0,  T].  The  infinite 
range  integral  can  be  used  within  this  time  frame  but  now  becomes  a  function  of  the 
signal  length  T,  as  well  as  the  frequency  /.  X{f,T)  is  the  finite  Fourier  transform  of 

*(*)•  j. 

X(/,r)=  /^*(t)e--'*'/‘*  (2) 

Jo 

The  aaalogne  signal  s(t)  can  be  sampled  at  N  discrete  time  intervals.  If  we  define  the 
sampling  interval  as  h  such  that  {N  -  l)h  =  T  and  h  is  chosen  to  produce  a  sufficiently 
high  cutoff  frequency,  the  above  finite  range  integral  can  be  expressed  as  a  summation. 
For  arbitrary  /, 

X{f,T)  =  h'£z„tTp  l-j2Tfnh]  (3) 

nsO 
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If  the  frequency  ipectnun  u  to  be  defined  »t  N  ducrete  freqnendei  (Le.  /  =  k/Nh), 
then  the  Diecrete  Fourier  Trantfarm  (DFT)  ii, 

X{fk,T)  =  h^‘*„exp[-i^l 

=  hXt  (4) 

where 


The  coiiq>utation  of  the  Fourier  components  Xk  at  the  discrete  frequency  values  is 
usually  carried  out  using  Fast  Fourier  Transform  (FFT)  techniques  which  are  discussed 
in  Section  3. 

Upon  calculation  of  the  Fourier  components,  it  is  possible  to  return  to  the  original 
input  signal  «„  by  performing  an  Inverse  Discrete  Fourier  Trsmsform  (IDFT).  The  IDFT 
equation  is: 

=  (6) 

Because  the  equations  for  OFT’s  and  IDFT’s  are  very  similar,  any  algorithm  used  to 
compute  a  DFT  can  also  be  used  to  calculate  its  inverse  by  singly  substituting  the 
discrete  frequencies  for  the  input  signal,  changing  the  sign  of  the  complex  argument 
from  negative  to  positive  and  finally  dividing  the  results  by  N. 

3  The  Fast  Fourier  Transform 

The  Fast  Fourier  lYansform  (FFT)  is  an  algorithm  that  enables  the  spectrum  of  an 
unsteady  signal  to  be  calculated  using  significantly  less  computational  effort  than  the 
standard  Fourier  series  procedure.  The  basis  of  the  FFT  is  to  split  an  N  point  transform 
into  two  N/2  point  transforms.  The  derivation  of  the  FFT  algorithm  from  the  DFT 
equation  (5)  is  s»  follows: 

Xk  =  -J-y- 

f»=0  ^  ^ 

=  53  I*, „ exp  +  »an+i  «P  j  | 

=  g  *jn  exp  +  «P  [-i^]  g  *>"+1  «*P 

This  can  be  re-written  as 

Xk  =  A%  +  Bf,Wf!,  (8) 

where 

=  51  *Jnexp|-y^^1 

n«0  ^  J 
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=  5^  *2n+l  «P  f-y— 

w^s  =  “p[-i^] 

As  an  exanqtle,  consider  a  case  where  the  input  signal  consists  of  eight  discrete  points. 
The  eight  point  transform  is  split  into  two  four  point  transforms  which  are  in  turn  split 
into  four  two  point  transforms.  The  equation  below  shows  this  splitting  process  with 
the  selected  data  points  identified  beneath  each  mdividual  summation 
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The  FFT  process  can  also  be  depicted  diagramatically  as  shown  for  the  eight  point 
transform  in  Figure  1  (Eef.  [2]).  The  order  in  which  the  raw  data  on  the  left  appears, 
is  determined  using  a  *bit  reversal’  procedure.  That  is,  the  address  of  the  n  th  point  is 
found  by  taking  the  mirror  image  of  the  binary  form  of  n.  For  example,  point  3  (Oil)  of 
the  8  point  transform  shown  is  relocated  to  position  6  (110).  The  dashed  lines  indicate 
that  the  quantity  be  multiplied  by  a  conqilex  coefficient  whereM  the  dotted  lines 
indicate  a  multiplication  by  unity.  The  junction  of  two  lines  at  the  dots  indicate  that 
the  two  quantities  should  be  added. 

When  calculating  the  Fourier  Transform  of  an  iV  point  signal  using  the  DFT  ap¬ 
proach,  it  involves  a  total  of  multiplications  each  followed  by  addition.  For  all  FFT 
applicsdions.  If  is  selected  to  be  a  power  of  2  (Le.  If  —  V  where  p  is  a  positive  integer), 
and  for  this  arrangement,  the  FFT  algorithm  reduces  the  number  of  computations  to 
W  log}  IV  (Ref.  |2]).  As  p  increases,  the  advantage  of  the  FFT  over  the  DFT  becomes 
apparent  as  indicated  in  Table  1. 


Table  1:  Conqtutational  requirements  for  Fourier  Transform  calculations 
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Figure  1:  TUuitration  of  the  FFT  proceii  for  an  JV  =  8  transform 


4  Spectral  Density 

The  Power  Spectral  Density  (PSD)  function  can  be  defined  as  the  rate  of  change  of 
the  mean  square  value  of  a  given  signal  vrith  respect  to  frequency  (Ref.  [3]).  The  mean 
square  value  of  an  infinitely  long  signal  z(t)  is  as  follows: 

<*> 

In  practice,  signal  analysis  is  confined  to  a  finite  length  time  history  x(t,  T)  where 

*(f,T)  =  *(t)  0<t<T 
=  0  otherwiie 

Therefore,  by  defining  z(t,  T)  between  ±oo,  the  mean  square  value  can  be  re-written  as 

=  lim  i  r  **(t,T)dt  (10) 

T-*oo  I  J-oo 

By  Tnalting  use  of  Parseval’s  theorem,  which  states  that  if  X(/)  is  the  Fourier  transform 
of  x{t)  then 

r  z[tydt=r  \x[f)fdf  (u) 

¥—00  V— oo 

and  since  X{f,T)  is  the  Fourier  transform  of  z{t,T),  the  mean  square  value  can  be 
defined  as  a  function  of  the  frequency  spectrum. 

=  2  lim  ^  r|X(/,r)|*<<r  (12) 

i  -•00  1  /o 
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The  PSD  function  <?z(/)  i*  the  rate  of  change  of  9^  with  respect  to  frequency,  therefore 

=  ^  =  [l-y(/.r)l*]  (13) 

As  ‘F’  denotes  the  expected  value  in  the  above  equation,  an  estimate  of  the  power 
spectral  density  function  can  be  expressed  as 

<?.(/)  =  I  W.r)!*  (14) 

with  the  tilde  (~)  indicating  that  the  function  is  an  estimate  only.  FVom  equation  (4) 
it  is  apparent  that  at  discrete  frequencies,  X(fi,T)  =  hXt.  The  power  spectral  density 
function  is  therefore  related  to  the  power  spectrum,  as  follows. 

ez(/ib)  =  ^|Xt|*  =  ^lX*l»  (15) 

For  a  second  signal  y(t),  another  PSD  estimate  can  be  calculated  using  the  above  equa¬ 
tions.  To  investigate  the  relationship  between  x(t)  and  y(t),  there  is  a  need  to  compute 
the  Cross  Spectral  Density  (CSD)  function.  The  definition  of  CSD  is  presented  below 
as  a  function  of  the  Fourier  components  of  each  signal. 

=  ~l^'(A,T)r(A,T)] 

=  (16) 

Here,  the  asterisk  (*)  indicates  the  complex  coitjugate. 

Discrete  frequency  spectra  estimates  are  subject  to  two  main  errors.  These  are 
discussed  in  Reference  [3]  and  are  summarised  as  follows. 

1.  Statistical  error  due  to  finite  data  length,  which  can  be  improved  by  averaging  the 
results  over  a  number  of  samples. 

2.  Bias  error  due  to  the  finite  bandwidth  filter  used  to  separate  the  various  frequency 
components  of  a  signal.  Bias  errors  are  large  when  the  PSD  changes  value  rapidly 
as  the  frequency  is  varied. 

5  Data  Manipulation 

5.1  Sampling  Considerations 

A  finite  length  signal  of  duration  T  should  be  sampled  at  JV  equally  spaced  time  intervals 
h.  The  parameters  AT  and  h  determine  the  range  and  resolution  of  the  final  frequency 
spectrum. 

5.1.1  Time  Interval,  h 

The  sampling  interval  h  is  selected  in  order  to  obtain  the  required  frequency  range.  If 
it  is  known  that  a  signal  contains  no  frequencies  greater  than  a  value  /c,  then  the  signal 


6 


can  be  completely  described  if  it  it  sampled  at  a  rate  greater  than  2/e.  Jc  is  known  as 
the  Nyquist  frequency  and  is  equal  to  half  of  the  sampling  frequency.  The  frequency 
spectrum  as  calculated  by  Fast  Fourier  Transform  techniques,  produces  a  mirror  inmge 
of  itself  about  tbis  point. 


(17) 


If  the  sampling  rate  is  selected  to  be  too  low,  superposition  of  the  high  and  low  frequency 
components  in  the  ori^nal  data  occurs.  Figure  2  shows  how  a  high  frequency  signal  is 
indistinguishable  from  a  lower  frequency  if  the  sampling  interval  is  too  large.  This  is 
called  ‘aliasing*. 


Figure  2:  The  aliasing  effect  in  the  time  domain 

For  any  frequency  /  within  the  range  0  <  /  </c.  the  higher  frequencies  which  are 
aliased  to  it  are 

(2/c±/),(4/c±/), . (2n/c±/) . 

Figure  3  illustrates  the  aliasing  effect  on  the  frequency  spectrum.  The  components 
that  lie  beyond  the  Nyquist  frequency  are  folded  back  about  this  point  to  corrupt  the 
true  frequency  spectrum. 

The  aliasing  phenomenon  can  be  avoided  by  either 

1 .  filtering  the  raw  analogue  signal  before  digitising  to  remove  frequencies  higher  than 
the  chosen  Nyquist  frequency  (i.e.  low  pass  filtering), 

2.  Selecting  h  small  enough  to  ensure  that  no  frequencies  with  significant  power  occur 
beyond  the  Nyquist  frequency. 
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Figure  3;  The  aliaiing  effect  in  the  frequency  domain 

5.1.2  Number  of  points,  N 

After  h  has  been  determined,  N  must  be  chosen  such  that  Nh  it  not  greater  than  the 
available  signal  length  T.  N  abo  controls  the  resolution  of  the  output  spectrum.  The 
number  of  discrete  frequency  points  within  the  range  (0,  fe)  i»  •Af/2  +  1.  JV  must  be 
selected  large  enough  to  enable  sudden  changes  (i.e.  spikes)  in  the  frequency  spectrum 
to  be  easily  identified,  and  to  avoid  this  information  becoming  lost  or  smeared  into 
neighbouring  frequencies.  The  number  of  points  must  also  be  selected  to  maintcun  a 
good  compromise  between  required  output  resolution  and  computational  effort. 

5.2  Centring 

Data  that  has  a  non-zero  mean  value  will  result  in  power  at  a  frequency  of  0  Hz.  If  this 
is  large  then  it  may  adversely  affect  the  results  at  other  frequencies.  The  time  history 
can  be  modified  to  centre  the  signal  on  a  zero  mean  value  to  remove  this  effect. 

*n  =  -  4  S  ** 


5.3  Windowing 

Because  it  is  impractical  to  compute  the  Fourier  transform  of  an  infinitely  long  signal,  it 
is  necessary  to  select  a  finite  length  signal  from  the  available  data.  This  is  accomplished 
by  multiplying  the  infinite  signal  by  a  boxcar  function,  u(t)  defined  as  follows. 

u(t)  =0  I  <  -T/2 

=  1  -T/2  <  t  <  T/2  (19) 

=  0  f  >  T/2 
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By  multiplying  the  infinite  signal  by  the  boxcar  function  u(t)  in  the  time  donuun,  the 
estimate  of  the  power  spectral  density  function  6x{f)  ^  been  affected  by  the  transform 
of  the  boxcar  function  If(/)  in  the  frequency  donuun. 

The  Fourier  transform  of  u(t)  is 

=  (20) 

and  is  called  a  tnndow  function  (Reference.  [1]). 

The  estimate  6x{f)  is  the  convolution  of  the  true  f?x(/)  with  the  window  function 
U{f).  The  main  lobe  of  the  window  function  spreads  the  power  at  discrete  frequencies 
over  the  sample  bandwidth  and  the  smaller  side  lobes  cause  leakage  to  occur  . 

To  reduce  the  effect  of  leakage,  it  is  necessary  to  modify  the  weighting  function  in  the 
time  domain  such  that  in  the  frequency  domain,  the  width  of  the  nuun  lobe  is  increased 
and  the  magnitude  of  the  side  lobes  is  decreased. 

One  common  technique  for  weighting  the  data  is  by  using  a  cosine  distributed  taper 
on  the  ends  of  the  signal.  This  ensures  that  there  are  no  discontinuities  of  amplitude 
and/or  slope  between  the  start  and  finish  of  each  sample  period  (Reference  [3]).  The 
Fourier  transform  is  based  on  the  assumption  that  the  data  is  periodic. 

If  the  total  number  of  digitised  points  in  the  signal  sample  is  N  and  the  number  of 
points  to  be  tapered  at  each  end  is  m,  then  the  weighting  function  Cn  for  n  =  0  to  iV  - 1 
is  given  by 


f  :n<m 

c„=M  :m<n<JV-m  (21) 

\l(l.co.(^))  :n>iV-m 

This  function  is  known  as  the  ‘Tukey’  window  (Ref.  [4]).  When  the  number  of  points 
at  each  end  is  equal  to  half  of  the  total  number  of  points,  the  Tukey  window  becomes 
commonly  known  as  the  ‘Hanning’  window.  Figure  4.  illustrates  the  difference  between 
the  Hanning  window  function  and  the  boxcar  window  function  in  the  frequency  domain. 

The  inclusion  of  a  taper  function  will  reduce  the  variance  of  the  tapered  data  with 
respect  to  the  original  data.  This  in  turn  reduces  the  magnitude  of  the  FFT  and  the 
subsequent  PSDs.  A  correction  sesde  factor  should  be  applied  to  the  resulting  PSD  if  a 
taper  function  has  been  used  in  order  to  compensate  for  the  reduced  variance.  Parseval’s 
theorem  suggests  that  this  factor  can  be  determined  by  integrating  the  square  of  the 
window  function  (i.e.  eq.22)  in  the  time  domain. 


•f  =  rm?dt 

J  — oo 

=  4/(4 -6^)  (22) 

Reference  [3]  indicates  that  for  a  full  Hanning  window  (i.e.  m/N  =  0.5),  the  above 
scale  factor  should  be  8/3  and  for  a  rectanguhr  window  (i.e.  m/N  =  0)  this  factor 
should  be  1. 
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Figure  4:  Comparison  of  Boxcar  and  Hanning  window  fimctiom. 

5.4  Smoothing  of  results 

With  real  data,  the  spectral  density  estimates  will  contiun  both  random  and  bias  errors. 
These  are  discussed  in  Reference  [1].  To  improve  the  accuracy  of  the  estimates,  it  is 
usual  to  average  the  results  using  one  of  the  following  techniques.  The  smoothed  result 
is  indicated  by  the  use  of  a  hat  ('')  which  replaces  the  tilde  (~)  used  to  indicate  a  raw 
estimate. 

The  time  history  (if  sampled  for  a  long  enough  period)  can  be  split  up  into  q  separate 
segments.  The  power  spectra  for  each  of  these  segments  should  be  calculated  indepen¬ 
dently  cmd  the  average  spectrum  determined  according  to  the  following  equation. 

G*  =  ^[<5M  +  <5*,J  +  --  +  6a^]  fc  =  0,l,2,...JV/2-l  (23) 

As  the  number  of  averages  increases,  the  results  swsynqitote  to  the  true  spectrum,  but 
there  is  a  practical  restriction  on  the  number  of  averages  that  can  be  used  due  to  the 
total  length  of  available  data.  In  order  to  increase  the  number  of  averages,  it  is  useful 
to  employ  an  ‘overlapping’  technique  where  the  tnuling  portion  of  data  in  the  current 
segment  can  be  used  as  the  initial  portion  for  the  next  segment.  An  overlap  of  more 
than  50  %  is  deemed  to  be  unnecessary,  as  each  new  data  block  then  becomes  too  similar 
to  the  last,  and  therefore  continued  avera^ng  would  yield  no  real  benefit. 

6  Auto-correlation  Functions 

The  auto-correlation  function  describes  the  relationship  between  a  signal  value  at  a 
particular  time  to  that  at  a  different  time.  For  a  sample  time  history  s((),  an  estimate 
of  the  auto-correlation  between  signal  values  at  times  t  and  t  +  r  can  be  made  by  taking 
the  product  of  the  two  signals  and  averaging  over  the  time  T.  This  process  is  represented 
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by  the  following  equation. 


=  z{t)*(t  +  T)dt  (24) 

Rx{t)  U  always  a  real  valued  even  function  (symmetric  about  r  =  0)  with  a  maximum 
at  r  =  0,  and  may  be  either  positive  or  negative  or  xero  at  all  other  values  of  r. 

Because  the  auto- correlation  function  can  also  be  defined  as  the  inverse  Fourier 
transform  of  the  power  spectrum,  it  can  be  efficiently  obt^ed  using  FFT  techniques. 

Rx(Tn)  =  (25) 

7  Cross-correlation  Functions 

In  a  similar  manner  to  the  auto- correlation  function,  the  cross-correlation  function  de¬ 
scribes  the  relationship  in  the  time  domain  between  two  different  signals.  For  a  pair  of 
signals  x(t)  and  y(t),  the  cross-correlation  function  is  defined  as: 

J  *(0»(f  +  r)dt  (26) 

ilzy(r)  is  a  real  function  that  can  be  positive  or  negative  or  sero  but  unlike  auto¬ 
correlation  functions,  it  does  not  necessarily  have  a  maximum  at  r  =  0  nor  does  it  have 
to  be  an  even  function. 

Again  by  using  FFT  procedures,  the  cross-correlation  function  can  be  calculated  as 
the  inverse  Fourier  transform  of  the  cross-spectrum. 

Averaging  of  a  number  of  ensembles  is  required  to  provide  estimates  with  satisfactory 
accuracy. 

Cross-correlation  plots  can  be  used  in  the  measurement  of  time  delays.  As  the  signal 
output  is  displaced  in  time  relative  to  the  input,  the  cross-correlation  plot  will  peak  at 
a  time  that  is  equal  to  the  time  difference.  It  should  be  noted  that  the  two  signals  must 
be  digitised  with  exactly  the  same  starting  time.  Any  misalignment  in  the  t  =  0  point 
will  result  in  a  phase  angle  error  in  the  cross  spectral  density  estimate,  which  will  in 
turn  cause  errors  in  the  cross-correlation  and  in  any  other  function  that  requires  CSV 
estimates  such  as  coherence. 

8  Coherence 

When  deaUng  with  two  signals  s(t)  and  y(t),  an  indication  of  the  relation  between  the 
signals  can  be  obt^ed  using  a  coherence  Action  7^(/)-  This  fuiction  varies  between 
0  and  1  over  the  specified  firequency  range.  When  is  equal  to  1,  the  signals  are 
siud  to  be  fiilly  coherent  (i.e.  there  is  a  linear  rdationship  between  them),  whereas  if 
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7|y  equali  0  the  (ignols  are  incoherent  or  totally  uncorrelated.  The  coherence  function 
relates  the  cross-spectrum  to  both  of  the  power  spectra. 


72,(A)  = 


5.(A)§,(A) 


(28) 


It  should  be  noted  that  the  above  equation  uses  the  hat  notation  indicating  that 
each  of  the  spectra  should  be  averaged  estimates.  The  use  of  raw  spectra  results  in  a 
coherence  function  that  is  everywhere  equal  to  unity,  since: 


-j  - .  ,  ^*»(A) _ 


(29) 


Smoothing  of  frequency  spectra  is  discussed  in  Section  5.4. 


9  Transfer  Functions 


If  we  consider  a  system  where  a  signal  p(t)  is  aresponse  to  an  inpnt  z(l),  the  relationship 
between  them  is  the  frequency  response  function  B{f).  For  a  single-input,  single-output 
system,  S(f)  (a  complex  function  of  frequency  which  contuns  both  gain  and  phase 
information)  can  be  best  calculated  using  the  averaged  cross  power  spectrum  and  the 
input  power  spectrum. 

S(/)=^  (») 

The  frequency  response  function  gun  is  the  the  modulus  from  the  above  equation. 


and  the  phase  angle  is  the  argument  of  S^y. 

?(/)  =  ar,(G„(/))  =  §^(/) 


(31) 


(32) 


10  Program  SPEC 

SPEC  (truncated  from  Spectral  Analysis)  is  written  in  standard  Fortran  77  except  for  a 
one  line  system  call  to  run  the  graphics  routine  GNUPLOT.  If  this  routine  is  unavailable 
or  unsuitable  for  the  computer  in  use,  then  the  line  should  be  commented  out  and  the 
program  recompiled. 

The  subroutines  used  by  SPEC  are  located  in  a  library  file  called  SPECLIB  which 
must  be  linked  to  SPEC  after  compilation.  The  library  includes  a  routine  to  remove 
the  mean  value  from  a  signal  (see  Section  5.2)  and  also  a  routine  to  apply  a  cosine 
wbidowing  function  to  the  raw  data  (see  Section  5.3).  The  most  important  routine 
within  this  library  is  the  FFT  routine  itself  which  was  obtained  from  Eeference  [5]. 
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10.1  Input  File  SPEC.IN 

SPEC  require!  the  creation  of  an  input  file.  This  file  contain!  aU  of  the  parameter! 
required  to  eet  up  the  calcuiation  and  ie  a  !iiiq>le  10  record  file  with  only  one  number 
on  each  record.  A  cample  of  !uch  an  input  file  ii  given  below.  The  comment!  within 
the  file  are  quite  acceptable  a!  SPEC  will  only  read  the  firet  number  on  each  line. 


1 

GinPLOT  plotting  to  bo  used  (  1  yea,  0  no) 

1 

Channel  nnsdser 

3 

Channel  number  (  0  it  only  one  channel  need) 

2000.0 

«« 

Stirling  rate  (Hz) 

300.0 

•  • 

Maximum  output  frequency 

30 

Ho.  of  averages 

10 

«* 

Exponent  for  no.  of  samples/average  (2*nl) 

60 

** 

Percentage  of  points  in  uindov  clip 

60 

Percentage  overlap 

0.09844 

Engineering  calibration  for  channel  1 

2.32616 

Engineering  calibration  for  channel  3 

The  purpose  of  each  record  is  as  follows. 

1.  Record  number  1  is  a  switch  to  turn  the  optional  GNUPLOT  plotting  routine  on 
or  off. 

2.  Record  numbers  2  and  3  select  the  data  files  to  be  used  in  the  analysis.  Here, 
channels  1  and  3  have  been  chosen  to  represent  signals  a  and  y  respectively  and 
these  datafiles  must  be  called  CHANOIJIAT  and  CHAN03.DAT.  If  record  number 
2  is  a  sero  then  the  analysis  will  be  for  a  single  channel  only  and  no  cross  signal 
calculations  will  be  attempted. 

3.  Record  number  4  is  the  rate  at  which  the  data  under  consideration  was  digitised 
from  the  analogue  signal. 

4.  Record  number  5  is  the  maximum  output  frequency  required  for  the  calculation 
which  cannot  be  greater  than  the  in«iTtiiim  Nyquist  cutoff  frequency  (half  of  the 
sampling  rate).  If  it  is  less  than  this  value  then  the  effective  sampling  rate  will  be 
reduced  by  reading  less  frequently  than  the  sampling  rate  to  satisfy  the  relation 
given  in  Section  6.1.1. 

5.  Record  number  6  is  the  number  of  segments  required  to  be  averaged.  The  pro¬ 
gram  will  a4iust  this  record  to  the  maximum  possible  if  there  is  not  enough  data 
available. 

6.  Record  number  7  is  the  base  2  logarithm  of  the  number  of  points  required  for  one 
segment.  For  the  above  example,  the  number  of  points  is  2*°  =  1024. 

7.  Record  number  8  is  the  percentage  of  the  total  number  of  points  required  at  each 
end  of  the  data  segment  for  the  cosine  window  clip.  For  the  above  exaiiq>le,  the 
number  of  points  is  0.5  x  2'°  =  612. 
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8.  Eecord  number  9  ii  the  percentage  overlap  of  raw  data  from  one  legment  to  the 
next. 

9.  Records  10  and  11  are  any  calibration  constants  that  need  to  be  applied  to  the  raw 
signal.  If  the  raw  signals  are  in  volts  and  they  need  to  be  converted  to  engineering 
units,  then  these  constants  should  be  specified  in  units/volt. 

10.2  Analysis  Format 

Assuming  that  the  input  file  SPEC  JN  is  available  and  is  in  the  correct  form,  the  spectral 
analysis  program  will  read  all  parameters  within  the  file  to  initialise  the  calculation.  A 
check  will  first  be  made  on  the  accessibility  of  the  requested  data  files  and  warning  flags 
will  be  output  to  the  screen  if  these  files  cannot  be  opened.  If  the  input  information  is 
correct,  a  summary  of  the  set-up  paranseters  is  provided.  This  allows  the  user  to  double 
check  that  the  analysis  has  been  specified  correctly. 

The  calculation  proceeds  with  the  reading  of  the  first  segment  of  data  from  either  one 
or  two  data  files.  The  data  segment  is  passed  to  a  centering  subroutine  (see  Section  5.2) 
and  then  the  cosine  windowing  function  is  applied  (Section  5.3).  The  data  is  adjusted 
according  to  the  calibration  constant  and  finally  the  resulting  real  number  is  converted 
to  a  complex  number  with  the  imaginary  component  equal  to  sero,  as  required  by  the 
Fast  Fourier  IVansform  algorithm.  The  results  from  the  FFT  routine  Xt,  Fj,  are  in  the 
form  of  complex  numbers  which  are  saved  in  an  array  while  the  next  data  segment  is 
read  and  processed  in  a  similar  manner.  The  following  parameters  are  calculated  and 
summed  over  the  total  number  of  data  segments. 

Xk,  n,  |X*|,  |F*|,  IXtl*.  |Ft|*,  XtF*,  \XkYk\ 

These  are  used  in  the  final  calculation  of  PSD  and  coherence,  etc.  After  the  requested 
number  of  averages  is  reached  or  the  total  available  data  has  been  read,  the  running 
totals  above  are  converted  into  averages  by  dividing  by  the  number  of  segments  used  in 
the  calculation.  A  summary  of  the  final  calculation  setup  is  then  provided  followed  by 
a  table  of  output  options. 

10.3  Output 

Selection  of  one  of  the  analysis  options,  will  create  a  data  file  into  which  the  results 
are  to  be  stored.  The  name  of  the  data  file  is  a  combination  of  the  type  of  analysis 
requested  (the  first  three  characters  in  the  file  name)  and  also  the  channel  numbers  used 
in  the  calculation.  For  example,  a  request  for  a  PSD  output  on  channel  1  would  have 
results  stored  in  file  PSDOl.OUT.  Similarly,  a  request  for  coherence  between  channels  1 
and  3  creates  a  file  COH01.03.OUT.  The  oupnt  files  possible  will  have  names  using  the 
following  codes: 

1.  DAT  -  Raw  data  from  the  last  segment  analysed. 

2.  FFT  -  The  average  Fourier  transform  over  the  time  history- 

3.  INV  -  The  inverse  Fourier  transform  of  the  average  FFT. 
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4.  PSD  -  The  average  Power  Spectral  Deniity.  Thii  option  also  calculates  a  variance 
estimate  by  integrating  the  PSD  function. 

6.  AUT  •  The  auto-ccBrelation  function  obt^ed  from  the  inverse  Fourier  transform 
of  the  average  PSD. 

6.  CSD  -  The  average  Cross  Spectral  Density. 

7.  CCO  •  The  cross-correlation  function  obtained  from  the  inverse  Fourier  transform 
of  the  average  CSD  (in  complex  form). 

8.  COH  -  The  coherence  function. 

9.  FRF  -  The  magnitude  of  the  frequency  response  function. 

10.  PHA  -  The  phase  angle  of  the  frequency  response  function. 

If  the  grsqthics  routine  GNUPLOT  has  been  made  available,  the  requested  output  will 
automatically  be  plotted  on  the  screen  for  instant  appraisal.  Further  analysis  can  be 
requested  until  the  EXIT  option  is  used.  All  data  files  created  by  SPEC  are  saved  for 
later  use  upon  exit  from  the  program. 

An  example  of  SPEC  output  is  presented  in  Appendix  A.  This  includes  a  comparison 
with  results  calculated  on  a  Wavetek  804A  Fourier  analyser. 

11  Conclusion 

The  development  of  a  PC  based  computer  program  (SPEC)  provides  a  capability  to 
perform  post  test  spectral  analysis  on  digitised  time-dependent  signals. 

The  mathematical  defuutions  have  been  included  to  ensure  that  the  basis  for  the 
calculations  used  in  the  computer  program  is  clearly  documented.  For  further  informa¬ 
tion  regarding  the  application  of  the  techniques,  reference  should  be  made  to  standard 
texts  such  as  Bef.  [1].  Some  general  guidelines  on  the  preparation  of  measured  data  are 
included. 
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Appendix  A  •  Example  SPEC  Output 

The  output  pretented  in  thii  appendix  hai  been  calculated  uiing  tignali  &am  an  ac¬ 
celerometer  and  an  untteady  pretture  transducer  obtained  during  testing  in  the  ABX 
9'  X  7'  Low  Speed  Wind  Tunnel.  Both  transducers  were  mounted  on  the  vertical  stabiliser 
of  a  l/9tb  scale  model  i^craft  which  was  subjected  to  aerodynamic  buffet. 

The  analysis  was  completed  under  the  following  conditions: 

1.  The  accelerometer  and  preuure  transducer  were  on  channels  1  and  2  respectively. 

2.  The  sampling  rate  was  2000  Hs. 

3.  Maximum  output  frequency  was  set  to  1000  Hs  (of  which  0  -  500  Hz  has  been 
plotted  in  the  following  figures). 

4.  30  block  averages  were  used  with  each  block  containing  1024  points  and  overlapping 
by  50  %. 

5.  The  data  was  unfiltered  at  the  time  of  recording  and  smoothed  using  a  iiill  Hanning 
window. 

A  typical  block  of  data  from  these  transducers  is  given  in  Figure  Al. 

Figure  A2  shows  the  correlation  functions  as  calculated  by  SPEC.  It  can  be  construed 
from  Parseval’s  theorem  that  the  value  of  autocorrelation  functions  at  r  =  0  should  be 
equal  to  the  variance  of  raw  signal.  Program  SPEC  provides  the  tuer  with  a  variance 
estimate  by  integrating  the  Power  Spectral  Density  function.  The  two  signals  chosen 
for  this  analysis  produced  variances  of  1376p^  and  lOSOOOPa^  which  compare  quite  well 
with  Figure  A2  at  r  =  0. 

The  Wavetek804A  Fourier  analyser  is  used  as  a  basis  for  comparison  with  SPEC 
output.  As  this  machine  does  not  have  a  built  in  facility  for  the  calculation  of  correlation 
functions,  a  direct  coiiq>aiison  with  a  hardware  based  Fourier  analyser  could  not  be 
made.  However,  for  the  parameters  presented  in  the  frequency  domain,  its  output  is 
useful 

For  the  frequency  based  fimctioiu  such  as  PSD,  coherence  etc.,  a  direct  comparison 
between  Wavetek  and  SPEC  results  is  graphically  presented  in  Figures  A3-A7.  This 
indicates  that  the  software-based  method  is  correct  in  determining  both  the  frequency 
information  within  a  random  signal  and  the  level  of  energy  at  each  discrete  frequency. 
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Tiiiie  sec. 


FIGURE  A2.  Coiielati<m  functions  calculated  using  program  SPEC 
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FIGURE  A3.  Power  Spectral  Density  of  Accelerometer  Signal 
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nCURE  A4.  Power  Spectral  Density  of  Pressure  Transducer  Signal 


FIGURE  A3.  Cross  Power  Spectral  Density  of  Above  Signals 
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FIGURE  A6.  Coherence  Function 
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FIGURE  A7.  Modulus  of  Transfer  Function 
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